function [Ireg, TBx, TBy, Fx, Fy, Mt] = affineregister_images(Imoving,Istatic,isaffine,iscomp,Options)
% This function register_images is the most easy way to register two
% images both affine and nonrigidly.
%
% Features:
% - It can be used with images from different type of scans or modalities.
% - It uses both a rigid transform and a nonrigid registration.
% - It uses multilevel refinement
% - It can be used with images of different sizes.
% - The function will automaticaly detect single modality or multiple
%   modalities, and choose the right registration method.
%
% [Ireg,Bx,By,Fx,Fy] = register_images(Imoving,Istatic,Options);
%
% Inputs,
%   Imoving : The image which will be registerd
%   Istatic : The image on which Imoving will be registered
%   Options : Registration options, see help below
%
% Outputs,
%   Ireg : The registered moving image
%   Bx, By : The backwards transformation fields of the pixels in
%       x and y direction seen from the  static image to the moving image.
%   Fx, Fy : The (approximated) forward transformation fields of the pixels in
%       x and y direction seen from the moving image to the static image.
%       (See the function backwards2forwards)
%   Mt : final affine transform matrix
% Options,
%   Options.SigmaFluid : The sigma of the gaussian smoothing kernel of the pixel
%                   velocity field / update field, this is a form of fluid
%                   regularization, (default 4)
%   Options.SigmaDiff : The sigma for smoothing the transformation field
%                   is not part of the orignal demon registration, this is
%                   a form of diffusion regularization, (default 1)
%   Options.Interpolation : Linear (default) or Cubic.
%   Options.Alpha : Constant which reduces the influence of edges (and noise)
%                   and limits the update speed (default 4).
%   Options.Similarity : Choose 'p' for single modality and 'm' for
%                   images of different modalities. (default autodetect)
%   Options.Registration: Rigid, Affine, NonRigid
%   Options.MaxRef : Maximum number of grid refinements steps.
%   Options.Verbose: Display Debug information 0,1 or 2
%
% Notes,
%   In case of Multiple Modalities affine registration is done with mutual
%   information. The non-rigid registration is done by first doing a
%   modality transformation (paints regions in image 1 with the intensity
%   pallette of those regions in image2 and visa versa), and than
%   using "normal" pixel based demon registration. See MutualTransform.m


% Function is written by D.Kroon University of Twente (March 2009)
% Modified by Tian Cao (March 2011)


%% add all needed function paths
try
    functionname='affineregister_images.m';
    functiondir=which(functionname);
    functiondir=functiondir(1:end-length(functionname));
    addpath([functiondir '/functions'])
    addpath([functiondir '/functions_affine'])
    addpath([functiondir '/functions_nonrigid'])
catch me
    disp(me.message);
end

% Disable warning
warning('off', 'MATLAB:maxNumCompThreads:Deprecated')

%% Process inputs
defaultoptions=struct('Similarity',[],'Registration','NonRigid','MaxRef',[],'Verbose',2,'SigmaFluid',4,'Alpha',4,'SigmaDiff',1,'Interpolation','Linear','maxiter',500);
if(~exist('Options','var')),
    Options=defaultoptions;
else
    tags = fieldnames(defaultoptions);
    for i=1:length(tags)
        if(~isfield(Options,tags{i})),  Options.(tags{i})=defaultoptions.(tags{i}); end
    end
    if(length(tags)~=length(fieldnames(Options))),
        warning('register_images:unknownoption','unknown options found');
    end
end

% Set parameters
MaxRef=Options.MaxRef;

% Start time measurement
if(Options.Verbose>0), tic; end

% Store the class of the inputs
Iclass=class(Imoving);

% Convert the inputs to double
Imoving=im2double(Imoving);
Istatic=im2double(Istatic);

% Resize the moving image to fit the static image
if(sum(size(Istatic)-size(Imoving))~=0)
    Imoving = imresize(Imoving,size(Istatic),'bicubic');
end

% Make smooth images for histogram and fast affine registration
ISmoving=imgaussian(Imoving,2.5,[10 10]);
ISstatic=imgaussian(Istatic,2.5,[10 10]);

%%Initialize Totoal Transfromation Field
%
%TBx = zeros(size(Imoving,1),size(Imoving,2));
%TBy = zeros(size(Imoving,1),size(Imoving,2));

%% Detect if the mutual information or pixel distance can be used as
% similarity measure. By comparing the histograms.
% if(isempty(Options.Similarity))
%     Hmoving= hist(ISmoving(:),60)./numel(Imoving);
%     Hstatic = hist(ISstatic(:),60)./numel(Istatic);
%     Hmoving(1)=0; Hstatic(1)=0;
%     if(sum(log(abs(Hmoving-Hstatic)+1))>0.3),
%         Options.Similarity='m';
%         if(Options.Verbose>0), disp('Multi Modalities, Mutual information is used'); drawnow; end
%     else
%         Options.Similarity='p';
%         if(Options.Verbose>0), disp('Same Modalities, Pixel Distance is used'); drawnow; end
%     end
% end
if(Options.Similarity(1)=='p'), type_affine='sd'; else type_affine='mi'; end

%% Register the moving image affine to the static image
% Affine register the smoothed images to get the registration parameters
if(strcmpi(Options.Registration(1),'R'))
    if(Options.Verbose>0), disp('Start Rigid registration'); drawnow; end
    % Parameter scaling of the Translation and Rotation
    scale=[1 1 1];
    % Set initial affine parameters
    x=[0 0 0];
elseif(strcmpi(Options.Registration(1),'A'))
    if(Options.Verbose>0), disp('Start Affine registration'); drawnow; end
    % Parameter scaling of the Translation, Rotation, Resize and Shear
    scale=[1 1 1 0.01 0.01 1e-4 1e-4];
    % Set initial affine parameters
    x=[0 0 0 100 100 0 0];
elseif(strcmpi(Options.Registration(1),'N'))
    if(Options.Verbose>0), disp('Start Affine part of Non-Rigid registration'); drawnow; end
    % Parameter scaling of the Translation, Rotation, Resize and Shear
    scale=[1 1 1 0.01 0.01 1e-4 1e-4];
    % Set initial affine parameters
    x=[0 0 0 100 100 0 0];
else
    warning('register_images:unknownoption','unknown registration method');
end

%pre affine transform
% isaffine = 0;
if isaffine>0
    for refine_itt=1:2
        if(refine_itt==2)
            ISmoving=Imoving; ISstatic=Istatic;
        end
        % Use struct because expanded optimset is part of the Optimization Toolbox.
        optim=struct('GradObj','off','GoalsExactAchieve',1,'Display','off','MaxIter',100,'MaxFunEvals',1000,'TolFun',1e-14,'DiffMinChange',1e-6);
        if(Options.Verbose>0), optim.Display='iter'; end
        x=fminlbfgs(@(x)affine_registration_error(x,scale,ISmoving,ISstatic,type_affine),x,optim);
    end
    
    % Scale the translation, resize and rotation parameters to the real values
    x=x.*scale;
    
    if(strcmpi(Options.Registration(1),'R'))
        % Make the rigid transformation matrix
        M=make_transformation_matrix(x(1:2),x(3));
    else
        % Make the affine transformation matrix
        M=make_transformation_matrix(x(1:2),x(3),x(4:5));
    end
end

%% Make center of the image transformation coordinates 0,0
[x,y]=ndgrid(0:(size(Imoving,1)-1),0:(size(Imoving,2)-1));
xd=x-(size(Imoving,1)/2); yd=y-(size(Imoving,2)/2);

% Initialize the modality transformed image variables
M_TF=[]; F_TF=[];


%% The nonrigid part of the registration
if(strcmpi(Options.Registration(1),'N'))
    
    % Demon registration parameters
    refinements=floor(log2(min(size(Imoving))/16));
    if(refinements>MaxRef), refinements=MaxRef; end
    parameters.sigma_diff=Options.SigmaDiff;
    
    % Non-rigid registration
    if(Options.Verbose>0), disp('Start non-rigid demon registration'); drawnow; end
    
    % Do every refinements step twice if modality transformation enabled
    if(Options.Similarity(1)=='m'),  loop=2; else loop=1; end
    
    %initialize ImovingTF, ImovingDF and Md
    %ImovingTF = Imoving;
    %ImovingDF = Imoving;
    %Tu = T;
    %Mt = id.
    Mt = [1 0 0; 0 1 0; 0 0 1];
    %initialize Mt with pre affine transform matrix
    if isaffine > 0
        Mt = M;
    end
    %Md = Mtu-Mt;
    %Q = id.
    Q = zeros(3,3);
    
    Bx = zeros(size(Imoving,1),size(Imoving,2));
    By = zeros(size(Imoving,1),size(Imoving,2));
    iteration = 0;
    %%
    %Igrid = imread('grid402.png');
    Igrid = imread('grid256.png');
    %Igrid = imread('grid.png');
    Ipre = Istatic;
    Idiff  = 100;
    
     while sum(abs(Q(1,1)-1)+abs(Q(1,2)-0)+abs(Q(2,1)-0)+abs(Q(2,2)-1)) > 0.001 %|| sum(abs(Md(1,3))+abs((Md(2,3))))> 1
        figure, hold on
        %while 1
        %for i=0:4
        %% Loop trough all refinements steps.
        %Bx = zeros(size(Bx,1),size(Bx,2));
        %By = zeros(size(By,1),size(By,2));
        % multiscale levels
        %         refinements = 0;
        for j=0:refinements
            %for j=0:1
            %             if refinements == 0
            %                 Bx = zeros(size(Bx,1),size(Bx,2));
            %                 By = zeros(size(By,1),size(By,2));
            %             end
            for l=1:loop
                % Set scaling parameters.resizepercentageentage
                resizepercentage=1/2^(refinements-j);
                if(resizepercentage>1), resizepercentage=1; end
                
                parameters.alpha=Options.Alpha*sqrt(resizepercentage);
                parameters.sigma_fluid=Options.SigmaFluid;
                
                if(Options.Verbose>0), disp(['Scaling resizepercentageentage : ' num2str(resizepercentage)]), end
                
                %% Incase of multiple modalities, transform both images to their
                % opposite modalities.
                %                 if(Options.Similarity(1)=='m')
                %                     if(Options.Verbose>0), disp('Start modality transformation'); drawnow; end
                %                     Bx_large=imresize(Bx,size(Imoving),'bicubic')*(size(Imoving,1)/size(Bx,1));
                %                     By_large=imresize(By,size(Imoving),'bicubic')*(size(Imoving,2)/size(By,2));
                %                     [Imoving_TF,Istatic_TF]=MutualTransform(Imoving,Istatic,15*sqrt(1/resizepercentage),4,Bx_large,By_large);
                %                     if(Options.Verbose>0), disp('Finished modality transformation'); drawnow; end
                %                 end
                %%
                sigma = 0.3/resizepercentage;
                % Set and resize the moving image and static image
                M = imresize(imgaussian(Imoving,sigma,[sigma*6 sigma*6]),resizepercentage,'bicubic');
                F = imresize(imgaussian(Istatic,sigma,[sigma*6 sigma*6]),resizepercentage,'bicubic');
                
                %% Resize the modality transformed images
                %                 if(Options.Similarity(1)=='m')
                %                     M_TF=imresize(imgaussian(Imoving_TF,sigma,[sigma*6 sigma*6]),resizepercentage,'bicubic');
                %                     F_TF=imresize(imgaussian(Istatic_TF,sigma,[sigma*6 sigma*6]),resizepercentage,'bicubic');
                %                 end
                
                %% Resize the transformation field to current image size
                Bx=imresize(Bx,size(M),'bicubic')*(size(M,1)/size(Bx,1));
                By=imresize(By,size(M),'bicubic')*(size(M,2)/size(By,2));
                
                % Put transformation fields in x and y direction in one variable
                B=zeros([size(M) 2]); B(:,:,1)=Bx; B(:,:,2)=By;
                % Store the dimensions of transformation field, and make a long vector from T
                sizes=size(B); B=B(:);
                
                % Parameters
                options.sigma_fluid=parameters.sigma_fluid;
                options.sigma_diff=parameters.sigma_diff;
                options.alpha=parameters.alpha;
                options.interpolation=Options.Interpolation;
                
                [B,Elist,Eulist] = demons_optimizer(M,F,Mt,B,sizes,iscomp,options);
                
                %% plot energy
                
                if resizepercentage > 0.99
                    if iteration == 0
                        Elistend = 0;
                        X = 1:5:size(Elist,2);
                        plot(X, Elist(1:5:end),'-ro' );
                        %plot(X, Elist(1:5:end),':b*' );
                        %plot(X, Eulist, '-g*');
                        %Elistend = size(Elist,2);
                    end
                    
                    if size(Elist,2)<=5
                        X = Elistend+1:Elistend+size(Elist,2);%+1;
                        plot(X, Elist,':b*' );
                        %plot(X, Eulist,'-g*')
                        Elistend = Elistend  + size(Elist,2);% + 1;
                    else
                        X = Elistend+1:5:Elistend+size(Elist,2);%+1;
                        plot(X, Elist(1:5:end),':b*' );
                        %plot(X, Eulist,'-g*')
                        Elistend = Elistend  + size(Elist,2);% + 1;
                    end
                    
                end
                
                % Reshape B from a vector to an x and y transformation field
                B=reshape(B,sizes);
                Bx=B(:,:,1); By=B(:,:,2);
                %%
            end
        end
        % do not do multiscale registrition
        if(refinements ~= 0)
            refinements = 0;
        end
        
        % Scale everything back if not already
        
        if(resizepercentage~=1)
            Bx=imresize(Bx,size(Imoving),'bicubic')*(size(Imoving,1)/size(Bx,1));
            By=imresize(By,size(Imoving),'bicubic')*(size(Imoving,2)/size(By,2));
        end
        
        %% Compute Affine part from the deformable transfrom
        
        DTX = Bx + xd;
        DTY = By + yd;
        
        %% Display results
        % affine transform displaement field in last iteration
        TFx = ((size(Imoving,1)/2) + Mt(1,1) * xd + Mt(1,2) *yd + Mt(1,3) * 1)-x;
        TFy = ((size(Imoving,2)/2) + Mt(2,1) * xd + Mt(2,2) *yd + Mt(2,3) * 1)-y;
        % composition of affine transform and deformable transform
        [TBx, TBy] = VFcomposition(TFx, TFy, Bx, By, size(Bx,2), size(Bx,1));
        
        % affine transform of moving image
        ImovingTF = displayImg(Imoving, TFx, TFy, 0); %title('Moving Image with Affine Transfromation');
        % deformable transform of the moving image after affine transform
        Itest1 = displayImg(ImovingTF,Bx,By,0); %title('Moving Image with Deformable Transformation');
        % composition transform of moving image
        Itestcomp = displayImg(Imoving,TBx,TBy,0); %title('Moving Image with Composition transform');
        
        % scale grid image
        Igrid = imresize(Igrid, [size(Istatic, 1) size(Istatic, 2)]);
        % affine transfrom of grid image
        Igridaffine = displayImg(Igrid, TFx, TFy, 0);
        % deformable transform of the grid image after affine transform
        Igridmov = displayImg(Igridaffine, Bx, By, 0);
        
        
        
        %% Demon's energy
        %
        %         iscomp = 1;
        if iscomp == 1
            Idiff_M1 = Itestcomp - Istatic;
        else
            Idiff_M1 = Itest1 - Istatic;
        end
        
        alpha = Options.Alpha;
        E_M1 =Idiff_M1.^2+alpha^2*Idiff_M1.^2.*(Bx.^2+By.^2);
        E1 = sum(E_M1(:));
        disp(['direct tranform energy = ', num2str(E1)]);
        
        %% SAD of fixed image and deformed moving image
        if iscomp == 1
            SAD1 = sum(sum(abs(Istatic-Itestcomp)));
            mse1 = sum(sum((Istatic-Itestcomp).^2))/(size(Istatic,1)*size(Istatic,2));
            ncc1 = NCC(Istatic, Itestcomp);
        else
            SAD1 = sum(sum(abs(Istatic-Itest1)));
            mse1 = sum(sum((Istatic-Itest1).^2))/(size(Istatic,1)*size(Istatic,2));
            ncc1 = NCC(Istatic, Itest1);
        end
        disp(['SAD1 = ', num2str(SAD1)]);
        disp(['MSE1 = ', num2str(mse1)]);
        disp(['NCC1 = ', num2str(ncc1)]);
        %         figure, imshow(Igridmov); title('Displacement field');
        %% plot vector field
        %figure, imshow(Imoving);
        %drawVectorField(Bx, By, x, y);
        
        %% get affine transformation
        % generate template
        mask = maskGen(Istatic, Istatic, 1);
        
        DTXmask = DTX.*mask;
        DTYmask = DTY.*mask;
        xdmask = xd.*mask;
        ydmask = yd.*mask;
        DTXmask(find(Istatic==0))=0;
        DTYmask(find(Istatic==0))=0;
        xdmask(find(Istatic==0))=0;
        ydmask(find(Istatic==0))=0;
        %[A, b] = GetAffineMatrix(DTX, DTY, xd, yd);
        %[A3, b3] = GetAffineMatrix2(DTX, DTY, xd, yd);
        [A, b] = GetAffineMatrix(DTXmask, DTYmask, xdmask, ydmask);
        
        %% solve similarity transform parameters
        DTmask = zeros(size(DTXmask,1)*size(DTXmask,2),2);
        xmask = zeros(size(xdmask,1)*size(ydmask,2),2);
        DTmask(:,1) = DTXmask(:);
        DTmask(:,2) = DTYmask(:);
        xmask(:,1) = xdmask(:);
        xmask(:,2) = ydmask(:);
        % get similarity transform parameters
        [pts, b1] = getAlignedPts( xmask, DTmask );
        
        %% get affine transform parameters
        %[A2, b2] = GetAffineMatrix2(DTXmask, DTYmask, xdmask, ydmask);
        
        Q = [A(1,1) A(1,2) b(1);
            A(2,1) A(2,2) b(2);
            0      0      1];
        Qinv = inv(Q);
        Mt = Mt*Qinv;
        % updated affine transform
        TFx = ((size(Imoving,1)/2) + Mt(1,1) * xd + Mt(1,2) *yd + Mt(1,3) * 1)-x;
        TFy = ((size(Imoving,2)/2) + Mt(2,1) * xd + Mt(2,2) *yd + Mt(2,3) * 1)-y;
        % affine transform
        Fx = ((size(Imoving,1)/2) + Qinv(1,1) * xd + Qinv(1,2) *yd + Qinv(1,3) * 1)-x;
        Fy = ((size(Imoving,2)/2) + Qinv(2,1) * xd + Qinv(2,2) *yd + Qinv(2,3) * 1)-y;
        % update the deformable displacement field
        Bx = ((size(Imoving,1)/2) + Q(1,1) * DTX + Q(1,2) * DTY + Q(1,3) * 1)-x;
        By = ((size(Imoving,2)/2) + Q(2,1) * DTX + Q(2,2) * DTY + Q(2,3) * 1)-y;
        % display deformable transform after extract affine transform part
        %figure, imshow(Bx,[]);
        %figure, imshow(By,[]);
        
        %% test
        %
        ImovingTF = displayImg(Imoving,TFx,TFy,0); %title('Updated Moving Image with Affine Transfromation');
        %drawDeformField(x+Fx, y+Fy, 'r');
        Itest2 = displayImg(ImovingTF, Bx, By, 0); %title('Affine Transformed Moving Image with Deformable Transformation');
        %ImovingTB = displayImg(Imoving,Bx,By,1); title('Updated Moving Image with Deformable Transfromation');
        %Itest2 = displayImg(ImovingTB, TFx, TFy, 1); title('Affine Transformed Moving Image with Deformable Transformation');
        
        % composition transform
        [TBx, TBy] = VFcomposition(TFx, TFy, Bx, By, size(Bx,2), size(Bx,1));
        [FBx, FBy] = VFcomposition(Fx, Fy, Bx, By, size(Bx,2), size(Bx,1));
        Itestcomp2 = displayImg(Imoving,TBx,TBy,0); %title('Updated Moving Image with Composited transformation');
        %drawDeformField(x+Fx, y+Fy, 'r');
        %figure,imshow(Istatic); title('static image'); colorbar
        %drawDeformField(x+FBx, y+FBy, 'r');
        
        %% demon's energy
        %Itest2 = Itest2.*Imask;
        if iscomp == 1
            Idiff_M2 = Itestcomp2 - Istatic;
        else
            Idiff_M2 = Itest2 - Istatic;
        end
        
        E_M2 =Idiff_M2.^2+alpha^2*Idiff_M2.^2.*(Bx.^2+By.^2);
        E2 = sum(E_M2(:));
        
        disp(['composited transform energy = ', num2str(E2)]);
        
        if iscomp == 1
            SAD2 = sum(sum(abs(Istatic-Itestcomp2)));
            mse2 = sum(sum((Istatic-Itestcomp2).^2))/(size(Istatic,1)*size(Istatic,2));
            ncc2 = NCC(Istatic, Itestcomp2);
        else
            SAD2 = sum(sum(abs(Istatic-Itest2)));
            mse2 = sum(sum((Istatic-Itest2).^2))/(size(Istatic,1)*size(Istatic,2));
            ncc2 = NCC(Istatic, Itest2);
        end
        
        disp(['SAD2 = ', num2str(SAD2)]);
        disp(['MSE2 = ', num2str(mse2)]);
        disp(['NCC2 = ', num2str(ncc2)]);
        %Igridmov = movepixels(Igrid, Bx, By, [], 3);
        
        %%
        iteration = iteration+1;
        if Idiff <= 1e-3
            break;
        else if sum(sum(abs(Q - eye(size(Q))))) < 1e-2
                break;
            else if iteration >= 100
                    break;
                end
            end
        end
        
        Idiff = sum(abs(Ipre(:)-Itest2(:)));
        Ipre = Itest2;
        %         Idiff = sum(abs(Ipre(:)-Itestcomp2(:)));
        %         Ipre = Itestcomp2;
        %end
        hold off
        %update Imoving
        %Imoving = ImovingTF;
        %end
        
        %TBy = ((size(Imoving,1)/2) + Mt(1,1) * Bx + Mt(1,2) * By + Mt(1,3) * 1)-x;
        %TBy = ((size(Imoving,2)/2) + Mt(2,1) * Bx + Mt(2,2) * By + Mt(2,3) * 1)-y;
        %TBx = TBx+Bx;
        %TBy = TBy+By;
        %% Transform the input image
        %Ireg=movepixels(Imoving,Bx,By,[], 3);
        
        %Affine transform
        TFx = ((size(Imoving,1)/2) + Mt(1,1) * xd + Mt(1,2) *yd + Mt(1,3) * 1)-x;
        TFy = ((size(Imoving,2)/2) + Mt(2,1) * xd + Mt(2,2) *yd + Mt(2,3) * 1)-y;
        %composited transform
        [TBx, TBy] = VFcomposition(TFx, TFy, Bx, By, size(TFx,2), size(TFx,1));
        
        % ImovingTF = movepixels(Imoving,TFx,TFy,[],3);
        %deformable transform
        % Ireg = movepixels(ImovingTF,Bx,By,[], 3);
        Ireg = movepixels(Imoving,TBx,TBy,[], 3);
        % Ireg = movepixels(Imoving,TFx+Bx,TFy+By,[], 3);
        % Ireg = Itest2;
        
        %debug info.
        disp(['SAD = ' num2str(sum(sum(abs(Istatic-Ireg))))]);
        
        
        %% debug info.
        %the composite transform is wrong
        %can not directly add two transform
        % TBx = TFx + Bx;
        % TBy = TFy + By;
        % [TBx, TBy] = VFcomposition(TFx, TFy, Bx, By, size(Bx,2), size(Bx,1));
        %for deformable field display
        % Bx = TBx;
        % By = TBy;
        %Itest = movepixels(Imoving, TBx, TBy, [], 3);
        %disp(['SAD = ' num2str(sum(sum(abs(Istatic-Itest))))]);
        
        %disp(['iteration = ' iteration]);
        disp([' Iteration = ' num2str(iteration)]);
        
        if ( nargout>3 )
            % Make the forward transformation fields from the backwards
            [Fx,Fy]=backwards2forwards(TBx,TBy);
        end
        
        % Set the class of output to input class
        if(strcmpi(Iclass,'uint8')), Ireg=im2uint8(Ireg); end
        if(strcmpi(Iclass,'uint16')), Ireg=im2uint16(Ireg); end
        if(strcmpi(Iclass,'uint32')), Ireg=im2uint32(Ireg); end
        if(strcmpi(Iclass,'int8')), Ireg=im2int8(Ireg); end
        if(strcmpi(Iclass,'int16')), Ireg=im2int16(Ireg); end
        if(strcmpi(Iclass,'int32')), Ireg=im2int32(Ireg); end
        if(strcmpi(Iclass,'single')), Ireg=im2single(Ireg); end
        
        % End time measurement
        if(Options.Verbose>0), toc, end
     end
end
end
%%
function [mask] = maskGen(Im1, Im2, scale)
mask = zeros(size(Im1,1),size(Im1, 2));
%         scale = 0.8;
wmask = size(mask, 2);
hmask = size(mask, 1);
offsetx = (1-scale)*wmask/2;
offsety = (1-scale)*hmask/2;
mask = Im1>0 & Im2>0;
tform = maketform('affine', [scale 0 0; 0 scale 0; offsetx offsety 1]);
mask = imtransform(mask, tform, 'bilinear',...
    'XData', [1 size(mask,2)],...
    'YData', [1 size(mask,1)]);

end


% function [FBx, FBy] = VFcomposition(AFx, AFy, BFx, BFy, width, height)
%
% [x,y] = ndgrid(0:height-1,0:width-1);
% xd = x - (height/2); yd = y - (width/2);
%
% AFTx = x + AFx;
% AFTy = y + AFy;
%
% BFTx = x + BFx;
% BFTy = y + BFy;
%
% FBx = zeros(height, width);
% FBy = zeros(height, width);
%
% xBas0 = floor(BFTx);
% yBas0 = floor(BFTy);
% xBas1 = xBas0 + 1;
% yBas1 = yBas0 + 1;
%
% xCom = BFTx - xBas0;
% yCom = BFTy - yBas0;
% perc0 = (1-xCom).*(1-yCom);
% perc1 = (1-xCom).*yCom;
% perc2 = xCom.*(1-yCom);
% perc3 = xCom.*yCom;
%
% %% limit indexes to boundaries
% check_xBas0 = (xBas0<0) | (xBas0>(height - 1));
% check_yBas0 = (yBas0<0) | (yBas0>(width - 1));
% xBas0(check_xBas0) = 0;
% yBas0(check_yBas0) = 0;
% check_xBas1 = (xBas1<0) | (xBas1>(height - 1));
% check_yBas1 = (yBas1<0) | (yBas1>(width - 1));
% xBas1(check_xBas1) = 0;
% yBas1(check_yBas1) = 0;
%
% Bx0 = AFx(1 + xBas0 + yBas0*height);
% Bx1 = AFx(1 + xBas0 + yBas1*height);
% Bx2 = AFx(1 + xBas1 + yBas0*height);
% Bx3 = AFx(1 + xBas1 + yBas1*height);
%
% By0 = AFy(1 + xBas0 + yBas0*height);
% By1 = AFy(1 + xBas0 + yBas1*height);
% By2 = AFy(1 + xBas1 + yBas0*height);
% By3 = AFy(1 + xBas1 + yBas1*height);
%
% %% set outside pixel to zero
% Bx0(check_xBas0 | check_yBas0) = 0;
% Bx1(check_xBas0 | check_yBas1) = 0;
% Bx2(check_xBas1 | check_yBas0) = 0;
% Bx3(check_xBas1 | check_yBas1) = 0;
%
% By0(check_xBas0 | check_yBas0) = 0;
% By1(check_xBas0 | check_yBas1) = 0;
% By2(check_xBas1 | check_yBas0) = 0;
% By3(check_xBas1 | check_yBas1) = 0;
%
% FBx = Bx0.*perc0 + Bx2.*perc2 + Bx1.*perc1 + Bx3.*perc3;
% FBx = reshape(FBx,width, height);
%
% FBy = By0.*perc0 + By1.*perc1 + By2.*perc2 + By3.*perc3;
% FBy = reshape(FBy,width, height);
%
% FBx = FBx + BFx;
% FBy = FBy + BFy;
%
% % FBx = FBx - x;
% % FBy = FBy - y;
%
% end

function Iout = displayImg(I, Bx, By, disp)
Iout = movepixels(I,Bx,By,[],3);
if disp == 1
    %         figure, imagesc(Iout);
    figure, imshow(Iout,[]); colorbar;
end
end

function drawDeformField(dx, dy, color)
hold on

contour( dx, 10, 'color', color );
contour( dy, 10, 'color', color );

hold off
end

function drawVectorField(Bx, By, x, y)
hold on
U = Bx;
V = By;
Xresize = imresize(x, 0.25);
Yresize = imresize(y, 0.25);
Uresize = imresize(U, 0.25);
Vresize = imresize(V, 0.25);
quiver(Xresize, Yresize, 15*Uresize, 15*Vresize);
%colormap hsv
hold off
end

function Ncc_val = NCC(Im1, Im2)
dIm1 = Im1 - mean(Im1(:));
dIm2 = Im2 - mean(Im2(:));

numerator = sum(sum(dIm1.*dIm2));
denominator = sqrt(sum(sum(dIm1.^2)))*sqrt(sum(sum(dIm2.^2)));

Ncc_val = numerator/denominator;

end